1 Introduction : la circulation à Paris

La ville de Paris rend publiquement accès à de nombreux jeux de données sur la circulation routière de la ville. La découverte du site web Paris Data (“Paris Data” visité le 26.12.2021) nous a motivé à utiliser cette ressource dans notre projet de prédiction. La possibilité de relier les différents capteurs dans une structure de graphe nous a intéressé. D’une part, cela motive des études du réseau entier : peut-on construire un algorithme global, qui apprend les interactions dans le réseau entier ? D’e l’autre côté’autre part, l’examination de rues individuelles mène à des questions sur la corrélation locale de la circulation. Un grand boulevard se comporte-t-il comme ses avenues voisines ? Enfin, dans tout cela, nous nous poserons la question à quel horizon temporel nous réussirons à prédire la présence de voitures dans les rues de Paris.

Avant de nous lancer, nous présentons d’abord le jeu de données en détail, sa richesse mais aussi les difficultés qu’il apporte : Sa taille et les valeurs manquantes. Ayant passé une grande partie du projet sur cet aspect, nous présentons nos stratégies pour résumer et compléter ce data set.

Cela fait, nous formulons plusieurs problèmes de prédiction, à différents horizons temporels et en utilisons différentes parties du data set. (TODO: spécifier à la fin pour avoir une petite table de matières ici)

2 Préparer le jeu de données

2.1 Les données et leurs attribus

Sur le site Paris Data, on trouve énormément de jeux de données issus des activités de la ville, notamment sur le comptage routier (“Données Du Comptage Routier” visité le 26.12.2021). Ces données sont collectées en continu grâce à des boucles électromagnétiques implantées dans la chaussée sur plus de 3000 tronçons de voies. Ces boucles et bornes de comptage mesurent principalement deux choses: le nombre de voitures qui passent et l’occupation du tronçon de route. Disponibles de 2014 jusqu’à aujourd’hui, ces deux données sont mesurées au pas horaire, résultant en environ \[ 8 \text{ ans } \times 365 \text{ jours } \times 24 \text{ heures } \times 3000 \text{ capteurs } = 2.1\times10^8 \text{ observations dans le data set.} \]

Chacune de ces lignes contient les informations suivantes, qui serons nos variables explicatives \(X_i\). Les noms des variables que nous utilisons dans notre analyse sont écrits en gras tout au long de cet exposé.

Un timestamp t_1h indiquant date et heure de la mesure.
Le nombre de voitures ayant passé le point de comptage pendant une heure, noté q et renommé par nbCar.
Le taux d’occupation, correspondant au temps de présence de véhicules sur la boucle en pourcentage d’une heure, noté k et renommé par rateCar.
La variable etat_barre indique si la route où se trouve la borne de comptage est “barée,” etat_barre=2, ou “ouverte,” 1. Nous la renommons state. Malheureusement, cette variable a souvent la valeur “inconnu,” 0, ou “invalide,” 4. La différence entre 3 et 4 est un mystère et par conséquence nous nous servons peu de cette variable.
L’identifiant iu_ac de chaque borne. Chaque ligne contient de informations sur les compteurs adjacents.
Le libellé (libelle) du tronçon de rue où est placé le capteur, par exemple “Champs Elysées partie B.”

Les données historiques de circulation sont enregistrées sur le site web Paris Data sous le format .txt. Chaque fichier couvre toutes les bornes pour une semaine. Pour commencer, il faut donc télécharger les fichiers de toutes les semaines et les convertir en dataframes. Ce faisant, nous effectuons également un tri par identifiant de borne iu_ac.

2.2 Simplification et aggrégation

Au vue de la quantité énorme de 200 millions d’observations, nous sommes obligés de réduire la taille du dataset. La première démarche consiste à garder les variables intéressantes: t_1h, nbCar, rateCar et iu_ac.

Ayant réduit la dimension de chaque observation, la deuxième étape est la réduction du nombre de capteurs de comptage. Pour ce faire, nous utilisons un modèle simplifié des rues de Paris. Etant donné que les capteurs couvrent les grandes axes routiers, dont le périphérique, notre objectif était d’agréger les données selon ces grands axes. Afin d’identifier où passent le plus de voitures, nous calculons les moyennes de nbCar pour chaque libellé, c’est-à-dire les moyennes à travers le temps et pour toutes les bornes associées à chaque libellé. En prenant les 200 libellés dont le nombre de voitures est le plus grand (hors périphérique), nous obtenons le graphique 2.1 qui représente très bien les axes routiers auxquels on s’attendait.
Ensuite, nous ignorons les libellés et regardons exclusivement ce graphique. De manière arbitraire, nous en déduisons le graphe de la figure 2.2, un schéma très simplifié de la circulation à Paris. Grâce à notre démarche d’utiliser les rues les plus fréquentées, cette abstraction nous permet d’avoir un modèle représentatif de la circulation parisienne.

Figure 2.1: Représentation en rouge des bornes d’observation associés aux 200 libellés les plus fréquentés (hors périphérique)

Figure 2.2: Graphe simplifié. Les points noirs indiquent les intersections entre nos arêtes.

Chaque arête dans le graphe 2.2 regroupe plusieurs capteurs dont nous agrégeons les données : après une collecte minutieuse et laborieuse des identifiants de toutes les bornes sur chacune des arêtes, nous associons pour chaque heure la moyenne de nbCar et rateCar des bornes de l’arête dont ils font partie. Cela a trois avantages:

Prendre la moyenne évite d’avoir des données manquantes, issues entre autres de capteurs dysfonctionnels sur une certaine période. Ce n’est pas toujours le cas, comme on le verra dans la section suivante.
S’il y a des capteurs où les données fluctuent beaucoup en raison de mauvais positionnement ou des perturbations indésirées, la moyenne réduit cette variance.
Le problème est réduit à 69 arêtes représentées par 69 dataframes.

Cette réduction des données résulte en la construction des jeux de données que l’on va manipuler et en la concrétisation du problème que l’on souhaite résoudre. Elle permet également d’obtenir des objets que nos ordinateurs sont capables de stocker efficacement. Il reste pourtant un problème: la présence de données manquantes.

2.3 Complétion de valeurs manquantes

En examinant des figures basiques des taux d’occupation, nous observons régulièrement une image comme dans le graphique TODO. Malgré l’agrégation le long des arêtes, il reste un nombre considérable de valeurs manquantes, NA, dans le dataset. Afin d’exécuter des tâches de prédiction, il nous faut cependant des données complètes. Une partie conséquente de notre travail a donc été de compléter les valeurs manquantes.

Figure 2.3: Taux d’occupation pour une journée ouvrière avec une valeur manquante

Figure 2.4: Visualisation de valeurs manquantes pour 30 arêtes (gris = NA)

En regardant le graphique 2.4, on considère qu’il y a deux types de “trous” dans le data set:

Des valeurs manquantes ponctuelles comme on a pu observer dans la figure 2.3 ou celles marquées en bleu dans le graphique ci-dessus. On peut deviner qu’il s’agit de simples erreurs dans la collection des données ou des travaux d maintenance par exemple.
Des trous plus grands, couvrant parfois plusieurs mois, comme ceux marqués en rouge dans le graphique 2.4. Il peut s’agir d’une absence étendue de plusieurs bornes, une clôture de la route en raison d’une construction, etc. Savoir si une route est effectivement inaccessible pendant de longues période serait une tâche de recherche en soi : comme mentionné plus tôt, la variable state ne contient presque aucune information. Nous remplissons alors les grands trous sans respecter d’éventuelles fermetures, ajoutant une autre simplification.

Le fait que les arêtes ont de différentes structures de valeurs manquantes peut aussi être observé en regardant la répartition du pourcentage de NA’s à travers les dataframes dans la figure 2.5.

Figure 2.5: Distribution des valeurs manquantes

Nous observons que les données de taux d’occupation sont généralement plus complètes : ceci est dû au fait qu’il y a plus de capteurs de ce type. Par ailleurs, nous avons vérifié qu’une valeur manquante de nbCar implique toujours un NA dans rateCar.

Pour remplir les trous, la première étape consiste à ajouter des variables explicatives supplémentaires. En fonction du type de trou, nous verrons qu’ils auront un impact important sur la complétion.

2.3.1 Variables supplémentaires : Proximité temporelle et locale

Pour l’instant, nous disposons des variables nbCar et rateCar pour chaque arête et chaque heure entre 2014 et 2020. Nous exploitons deux propriétés structurelles des données pour obtenir de l’information supplémentaire:

D’abord, il s’agit de séries temporelles, donc nous rajoutons les valeurs historiques de nbCar et rateCar décalées d’une heure, d’un jour et d’une semaine (nbCarLaggedHour, nbCarLaggedDay, et nbCarLaggedWeek) comme variables explicatives. L’espoir étant que la circulation reste identique à travers ces cycles. Avec l’objectif spécifique de compléter les données manquantes, il est aussi pertinent d’ajouter ces valeurs du futur (rateCarFuturHour etc.), étant donné que le modèle d’imputation a pour but d’interpoler au lieu d’extrapoler. On fait l’hypothèse que toutes ces variables décalées seront efficaces pour remplir les trous locaux : une heure ou une journée manquante devrait facilement être inférée en utilisant des données qui sont proches dans le temps.

Pour les valeurs manquantes au milieu des grands trous, nous n’avons pas accès aux heures et jours d’avant car elles sont aussi manquantes. Cependant, grâce à l’interconnexion de nos données à travers le graphe présenté plus haut, les mesures de circulation d’une arête peuvent bien être expliquées en fonction de celle de ces voisines. Dans un effort manuel, nous avons pour chaque arête \(A\) collectionné des voisins \(V_i\), dans le sens que les voitures dans \(A\) passent ensuite par l’un des \(V_i\) ou bien l’inverse. Nous faisons ici encore des choix arbitraires de qui est voisin de qui pour deux raisons. Premièrement, inférer statistiquement quelle arête influence quelle autre nécessiterait des données déjà complètes pour une régression par exemple. Deuxièmement, utiliser toutes les autres 68 arêtes au lieu de juste 2 à 6 voisins augmenterait trop le temps d’exécution de la complétion.
Donnons un exemple illustrant le choix de voisins ainsi que la dénomination des nouvelles variables. Parmi les voitures qui entrent dans Paris par l’arête “pont amont - pont austerlitz,” il y a une partie considérable qui continue tout droit sur l’arête “pont austerlitz - chatelet.” Par conséquent, nous ajoutons la variable rateCar_pontausterlitzTOchatelet au dataframe de “pont amont - pont austerlitz.”

Jusqu’ici, nous n’avons utilisé que des informations déjà présentes dans le jeu de données. Dans la suite, nous ajoutons des variables extérieures qui pourraient également expliquer la circulation routière.

2.3.2 Variables supplémentaires : Facteurs externes

Variables temporelles

Le traffic routier étant relié à l’activité humaine, nous avons ajouté de nombreuses variables temporelles (notamment grâce à la librairie lubridate).

year, month, day, hour résultent d’une simple décomposition des timestamps t_1h,
time est un simple compteur de journées, commençant à \(1\) le 1er janvier 2014,
toy (time of year), il s’agit d’un numéro entre 0 et 1 indiquant la position de l’observation dans l’année en cours,
weekdays le numéro du jour de la semaine et weekendsIndicator l’indicatrice si le jour est un jour du weekend,
winterHolidaysIndicator et summerHolidaysIndicator les indicatrices des vacances d’hiver et d’été définies à partir de (“Les Archives Du Calendrier Scolaire” visité le 26.12.2021),
bankHolidaysIndicator l’indicatrice des jours fériés définie à partir de (“Jours Fériés” visité le 26.12.2021).

Index de la situation sanitaire en rapport avec le Covid-19

Nous avons également ajouté, covidIndex, un index allant de 0 à 100 représentant les restrictions du gouvernement sur la situation sanitaire en rapport avec le Covid-19. Il est calculé à partir de nombreux indicateurs et est fourni par l’université d’Oxford (“Covid Index” visité le 26.12.2021). Cette variable permettra éventuellement une étude de l’année 2020 qui voit une circulation perturbé. Pour cela, il faudra utiliser des modèles qui s’adaptent très vite à l’influence de nouvelles variables.

Météo

A partir de données de l’Organisation Météorologique Mondiale (“Observations Météorologiques Historiques En France” visité le 26.12.2021), nous avons extrait 2 variables météorologiques : temperature la température en Kelvin et precipitation les précipitations dans les 3 dernières heures en mm. Ces données ont été mesurées à Athis-Mons en Essonne.

On dispose d’un relevé tous les 3 heures environ donc on procède à une interpolation pour compléter les données. Etant donné que les mesures sont uniformément réparties, on utilise une interpolation linéaire basique à l’aide de la fonction na.approx de la librairie zoo.

2.4 Complétion via Random Forest

Avec un dataframe enrichi à notre disposition pour chacune des arêtes du graphe de circulation, nous pouvons passer au remplissage des valeurs manquantes. Nous sommes pourtant restreint dans notre recherche d’un algorithme d’imputation car plusieurs variables explicatives contiennent également des NA. D’une part, nous voulons remplir deux variables en même temps (nbCar et rateCar). D’autre part, les voisins n’ont pas moins de valeurs manquantes, laissant de grands trous dans le dataframe.
Heureusement, il y a un modèle qui peut être entraîné en dépit de valeurs manquantes parmi les variables : les forêts aléatoires constituées d’arbres CART. Une règle de décision dans un tel arbre peut être ignorée si la valeur de la variable nécessaire n’est pas renseignée. Le package miceRanger profite de ce fait en utilisant un algorithme dit d’imputation multiple : en commençant par la variable \(V_1\), une forêt aléatoire est entraînée sur les lignes où \(V_1\) n’est pas manquant. Puis, les autres lignes sont “prédites” par la forêt. Avec le nouveau dataset moins vide, le processus est reproduit et ainsi de suite (Van Buuren 2018). Nous nous arêtons pourtant à la deuxième itération vu qu’il n’y a pas d’intérêt à compléter les données des voisins.

Comme nous répétons le processus de complétion 69 fois, le temps d’exécution est de quelques heures sur nos ordinateurs. Nous n’avons donc pas le luxe d’optimiser des hyperparamètres comme nous ferons plus tard pour la tâche de prédiction. Le choix de 100 arbres par forêt et 7 variables considérées à chaque split semble un bon compromis. TODO Guillaume: Voulons-nous mettre du code ici?
Au délà des capacités de complétion, les forêts aléatoires permettent de calculer des scores d’importance pour chaque variable explicatives. Ces scores peuvent aider à identifier lesquelles des variables ont le plus contribué à l’estimation de la cible (le score est calculé comme la réduction de variance à chaque split). Pour une arête du graphe, jussieu - saint michel, nous illustrons les variables qui ont le plus contribué à l’imputation :

Figure 2.6: Variables avec le score d’importance le plus haut pour rateCar de l’arête ‘’jussieu - saint michel’’

Figure 2.7: Variables avec le score d’importance le plus haut pour nbCar de l’arête ‘’jussieu - saint michel’’

Sur l’arête dont traitent les deux plots, on observe que rateCar et nbCar exhibent deux comportement assez différents. Pour le taux d’occupation, ce sont presque exclusivement les valeurs temporellement décalées qui contribuent à l’imputaion. Comme il y seulement 0.8% de valeurs à compléter, ceci confirme notre hypothèse que les petits trous d’une ou plusieurs heures sont très bien complétés par les valeurs des heures juste avant ou après. Pour le nombre de voitures, où la proportion de NA était plus importante, les deux variables ayant le plus réduit les variances aux splits sont des variables nbCar issues d’arêtes voisines. Les trous importants en nbCar entre Jussieu et Saint Michel semblent être le mieux expliqué par ce qui se passe autour. Cela confirme notre choix de variables pour combler les 2 types de trous. Un phénomène similaire peut s’observer dans les autres arêtes.(TODO Guillaume: Quoi d’autre écrire sur ceci ? Y a des scores intéressants dans RF (genre error Out of Bag) que j’aurais dû retenir de l’imputation ?)

Après l’imputation arête par arête, nous pouvons enfin compléter toutes les données grâce au fait que les variables reliées aux voisins ont été imputées dans les dataframes qui correspondent à ces voisins. Il est important de noter ici que nous allons dans la suite couper le jeu de données en deux (train - test) et que cette imputation sera effectuée individuellement sur chacune des deux parties. Faire ceci deux fois séparément est nécessaire pour préserver l’indépendance et permet ensuite d’estimer les erreurs de nos modèles.
Dans la suite, nous supposons que les données sont complètes. Il faut pourtant remarquer que ces forêts aléatoires font indirectement partie des modèles de prédiction que nous allons utiliser. De plus, nous gardons à l’esprit que nous avons introduit un biais dans nos données d’apprentissage comme de test.

3 Modèles de prévision

3.1 Définition des problèmes

Avec un jeu de données complet, nous pouvons enfin utiliser des algorithmes de Machine Learning pour faire de la prévision. Vu la complexité du jeu de données, il n’est même pas évident quelles questions se poser. Il y a deux dimensions principales du problème que nous allons examiner dans la suite, en utilisant de différents modèles d’apprentissage:

La première dimension du problème est temporelle : Quel moment dans le futur voulons-nous prédire ? Autrement dit, jusqu’à quel horizon temporel est-ce que nous nous permettons d’utiliser le passé pour prédire le présent ?
1. Le cas une heure correspond à un problème du genre “Google Maps” : si nous prenons la voiture maintenant, quelle sera la circulation dans les prochaines 60 minutes ?
2. Etant donnée la circulation du jour précédent, la tâche devient intéressante pour les régulateurs des feux de circulation ou pour la police qui pourront ainsi éviter des bouchons.
3. Si nous voulons prédire une semaine voire un mois en avance, le modèle devra comprendre davantage le comportement habituelle de la circulation routière. Ce sera plutôt un pari basé sur des fluctuations habituelles
La deuxième dimension est reliée au réseau routier : Pouvons-nous profiter de la structure du graphe pour améliorer la prédiction ?
1. La question de base sera si on réussit à battre des benchmarks sur une seule arête.
2. Puis nous pourrons reutiliser les “voisins” introduits pendant l’imputation.
3. Finalement, une question est si nous réusissons à trouver un seul modèle décrivant toutes les arêtes en même temps, une sorte de prédiction multivariée.

Avant de nous lancer, nous précisons comment nous évaluons nos modèles et quels sont les benchmarks, les modèles naïfs, que nous voulons battre.

3.1.1 Evaluation des modèles

Train - test - split

On découpe notre jeu de données en 2 parties de proportion 2/3 et 1/3 : une partie apprentissage de 2014 à 2017 et une partie test de 2018 à 2019. Faire le découpage de cette manière assure qu’on est évalué sur les prédictions du futur plutôt que du passé. L’année 2020 est omise à cause du Covid, mais elle pourrait faire partie d’un autre projet sur le jeu de données en relation avec la pandemie !

Mesurer la performance

Pour évaluer les performances des modèles, nous considérerons le Root Mean Square Error (RMSE) qui représente la racine carrée du deuxième moment d’échantillonnage des différences entre les valeurs prédites et les valeurs observées. Une autre métrique souvent utilisée est le Mean Absolute Percentage Error (MAPE) qui n’est pas applicable à notre situation car un nombre important des valeurs de rateCar sont nulles et on ne peut pas diviser par elles.

Modèles naïfs à battre

Afin de comparer nos modèles, on construit 3 modèles témoins. Ces modèles sont dit naïfs car extrêmement simpliste :

naiveModel1 prévoit le taux d’occupation d’une heure \(t \in \{ 0, \dots, 23 \}\) d’un jour \(j \in \{ 1, \dots, 7\}\) d’un mois \(m \in \{ 1, \dots, 12 \}\) d’une année \(a \in \{ 2018, 2019 \}\) en moyennant les taux d’occupation de l’heure \(t\) du jour \(j\) du mois \(m\) pour toutes les années \(a \in \{ 2014, \dots, 2017 \}\).
naiveModel2 fait de même en calculant la médiane et non la moyenne.
naiveModel3 est simplement le taux d’occupation à l’heure précédente rateCarLaggedDay.

On calcule la moyenne des erreurs RMSE pour les 69 arêtes :

Table 3.1: Performances des modèles naïfs témoins
naiveModel1	naiveModel2	naiveModel3
5.844	5.899	3.955

Les 2 premiers modèles naïfs ont des scores similaires. Le troisième modèle naïf, qui consiste à prendre la variable rateCar lagged, est meilleur. On interprète cette bonne performance par la très forte corrélation du taux de circulation d’une heure avec la suivante. Dans la suite, il sera intéressant de comparer ce modèle naïf avec des modèles plus complexes pour justifier leur utilisation.

Recherche de paramètres

Afin d’optimiser les performances de nos différents modèles, nous effectuons de la recherche de paramètres. On utilise les 2 méthodes suivantes :

Validation croisée sur 16 blocs. On divise ainsi nos 4 années de données d’apprentissage en blocs de 3 mois, ce qui englobe les différentes saisonnalités de nos données. On a codé cette méthode à la main dans la fonction ….
Validation croisée progressive, où on fixe une fenêtre initiale d’apprentissage de 2 ans, que l’on incrémente dans l’ordre chronologique de 2 mois à chaque itération, pour entraîner le modèle et mesurer sa performance à l’horizon 1 (c’est-à-dire prédire le taux d’occupation de l’heure suivante). Pour utiliser cette méthode, on utilise le package caret avec le paramètre timeSlice.

Validation croisée progressive [@CVprogressive]

Figure 3.1: Validation croisée progressive (“Progressive Cross-Validation” visité le 22.02.2022)

La caractéristique des données de type série temporelle est qu’elles sont ordonnées et possèdent une dépendance temporelle. C’est pourquoi il semble pertinent de ne pas mélanger les observations avant de découper le jeu de données en données apprentissage et données test. De plus, il semble aussi important de ne pas utiliser des données futures pour prédire le passé. Dans la validation croisée, il est donc naturel de retirer les blocs d’apprentissage ultérieurs au bloc test. C’est la validation croisée progressive. Cependant, en pratique, on souhaite entraîner le modèle une unique fois sans le mettre à jour tous les 2 mois. Dans ce cas, la validation croisée est plus pertinente.

Lors de nos recherches de paramètres, nous avons observé que les 2 méthodes sélectionnent des paramètres proches. Finalement, on se base sur la validation croisée en 16 blocs pour la recherche de paramètres.

Pour permettre la reproductibilité des résultats, on fixe la graine du générateur aléatoire. De plus, on effectue la recherche de paramètre pour une seule arête sur les 69 car c’est un processus long et on fait l’hypothèse que le paramètre optimal pour la prévision sur une arête est proche de celui de toutes les autres arêtes.

3.2 Google Maps : Prédire la prochaine heure

Commençons par la tâche qui est censée être la plus simple : Prédire la circulation dans les prochaines 60 minutes.

3.2.1 GAM pour une arête

Pour l’instant, nous nous restreignons à une arête isolée dans le sens que nous ne prenons pas en compte la circulation sur ces voisins. Notre objectif initial était de faire une simple régression linéaire, mais il y a quelques relations pas linéaires du tout. Si on essaie par exemple, d’écrire le taux d’occupation d’une avenue en fonction de l’heure, le graphique suivant montre bien que les heures de pointe seraient mieux décrites par un polynôme que pour une droite :

Figure 3.2: Une journée typique sur l’arête

Pas conséquent nous travaillons plutôt sur des GAM que sur de simples régression. Constituer un modèle additif et en sélectionner le meilleur est une tâche sans fin. C’est pourquoi nous adapterons une approche greedy pour la sélection de modèles où nous rajoutons une variable explicative à la fois. Idéalement ce sera celle qui apporte le plus de puissance prédictive supplémentaire (forward search).

Pour prédire une heure dans le futur, il est incontestable que l’heure présente contient le plus d’informations. La variable rateCar_LaggedHour sera donc la première à considérer. L’effet de l’occupation d’une heure à l’autre n’est pourtant pas homogène à travers les heures de la journée. En plus, comme nous voyons dans le plot supérieur gauche, il y a peut-être une relation linéaire entre les variables, qui a pourtant une grand variance :

Les autres trois plots que nous venons de voir représentent la même comparaison de rateCar et de rateCar_LaggedHours, mais sur différentes heures. On a bien l’impression que des régressions linéaires pourront nous servir, vu la variance qui est visiblement réduite. Comparer le nuage de points à 8 heures avec ceux de 12 et de 13 heures nous donne envie de regrouper certaines heure. Ce phénomène était déjà observable dans le graphique 3.2 : dans l’heure de pointe matinale (morning rushhour), la circulation augmente beaucoup entre 7 et 8 heures tandis qu’entre 12 et 13 heures elle reste plutôt constante.

Afin de regrouper les heures de la journées en “périodes,” nous procédons de manière heuristique : pour chaque heure de la journée, nous effectuons une régression linéaire et comparons les coefficients (intercept et slope). Ceci donne le tableau suivant.

##  hour intercept slope
##     0    -0.148 1.050
##     1     0.038 0.749
##     2    -0.218 0.736
##     3     0.039 0.651
##     4     0.007 0.645
##     5     0.163 0.578
##     6     0.414 0.575
##     7     1.155 0.568
##     8    -0.630 2.993
##     9    -1.185 2.298
##    10     1.691 0.846
##    11     2.253 0.888
##    12     1.688 0.949
##    13     2.278 0.879
##    14     0.756 0.808
##    15     1.778 0.842
##    16     2.442 1.046
##    17     1.527 1.013
##    18     1.230 0.916
##    19     1.529 0.810
##    20     2.020 0.795
##    21     0.709 0.685
##    22     2.252 0.274
##    23     1.675 0.438

Les résultats sont congruents avec le graphique 3.2, par exemple l’explosion de slope à 8 heures où commence l’heure de pointe matinale. En nous basant sur ce tableau, nous regroupons de manière heuristique les heures en différents groupes: “night” de 0 à 7 heures, “8heures” à huit heures, “9heures” similairement, “noon” de 10 à 15 heures, “afternoon” de 16 à 17 heures, “evening_rush” de 18 à 20 heures et “evening” de 21 à 23 heures. Ainsi, chaque ligne du dataframe reçoit le nouveau attribut hour_groups. Avec ces groupes, nous pouvons établir notre premier modèle mod1, une régression linéaire avec un “mixed effect” hour_groups * rateCar_LaggedHour :

cluster_hours = function(h){
  if (h %in% c(0,1,2,3,4,5,6,7)){return("night")}
  if (h %in% c(8)){return("8heures")}
  if (h %in% c(9)){return("9heures")}
  if (h %in% c(10,11,12,13,14,15)){return("noon")}
  if (h %in% c(16,17)){return("afternoon")}
  if (h %in% c(18,19,20)){return("evening_rush")}
  if (h %in% c(21,22,23)){return("evening")}
  else{return("asfd")}
}
df_train$hour_groups = sapply(as.factor(df_train$hour), cluster_hours)
df_test$hour_groups = sapply(as.factor(df_test$hour), cluster_hours)

# for (group in unique(df_train$hour_groups)){
#   indices = df_train$hour_groups == group
#   plot(df_train[indices,]$rateCar_LaggedHour, df_train[indices,]$rateCar,
#        main = group, xlab = "rateCar_laggedHour", ylab = "rateCar")
# }

form1 = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour
mod1 = lm(data = df_train, formula = form1)
summary(mod1)

## 
## Call:
## lm(formula = form1, data = df_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.480  -0.836  -0.301   0.749  51.016 
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                0.846441   0.024047   35.20   <2e-16
## hour_groups8heures:rateCar_LaggedHour      2.201884   0.038810   56.73   <2e-16
## hour_groups9heures:rateCar_LaggedHour      1.935035   0.013804  140.18   <2e-16
## hour_groupsafternoon:rateCar_LaggedHour    1.097172   0.003656  300.10   <2e-16
## hour_groupsevening:rateCar_LaggedHour      0.572480   0.003810  150.26   <2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.881992   0.002559  344.71   <2e-16
## hour_groupsnight:rateCar_LaggedHour        0.563815   0.011286   49.96   <2e-16
## hour_groupsnoon:rateCar_LaggedHour         0.928103   0.002810  330.23   <2e-16
##                                               
## (Intercept)                                ***
## hour_groups8heures:rateCar_LaggedHour      ***
## hour_groups9heures:rateCar_LaggedHour      ***
## hour_groupsafternoon:rateCar_LaggedHour    ***
## hour_groupsevening:rateCar_LaggedHour      ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour        ***
## hour_groupsnoon:rateCar_LaggedHour         ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.421 on 35052 degrees of freedom
## Multiple R-squared:  0.8753, Adjusted R-squared:  0.8753 
## F-statistic: 3.516e+04 on 7 and 35052 DF,  p-value: < 2.2e-16

Toutes les variables ont une contribution significative d’après les tests de student. Calculons le score de validation croisée pour ce modèle simple pour obtenir

CV_gam = function(formula, df_train, no_folds){
  K = no_folds 
  data = df_train
  folds <- cut(seq(1,nrow(data)),breaks=K,labels=FALSE)
  fold_scores = c()
  start_time = Sys.time()
  for(i in 1:K){
    testIndexes <- which(folds==i,arr.ind=TRUE)
    testData <- data[testIndexes, ]
    trainData <- data[-testIndexes, ]
    g = gam(data=testData, formula=formula)
    Y_test = testData$rateCar
    Y_predict = predict(g, newdata=subset(testData, select = -c(rateCar)))
    fold_scores = c(fold_scores, rmse(Y_test, Y_predict))
  }
  exec_time = Sys.time()-start_time
  score = mean(fold_scores)
  ret = data.frame(score = score, no_folds = no_folds, exec_time = exec_time, 
                   formula = paste(as.character(formula)[c(2,1,3)], collapse = " "))
  return(ret)
}

CV_gam(form1, df_train, no_folds = 8)

##     score no_folds      exec_time
## 1 2.39375        8 0.8233941 secs
##                                                                         formula
## 1 rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour

Afin de choisir la prochaine variable à rajouter pour expliquer rateCar, nous les explorons une par une :

hour :

form1_hour = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + 
  s(hour, bs="cc", by=weekendsIndicator)
# mod1_hour = gam(form1_hour, data=df_train)
# summary(mod1_hour)
CV_gam(form1_hour, df_train, no_folds = 8)

##     score no_folds    exec_time
## 1 2.29125        8 2.03911 secs
##                                                                                                                      formula
## 1 rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(hour, bs = "cc", by = weekendsIndicator)

weekdays :

form1_hour = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(hour, by=weekendsIndicator)
# mod1_hour = gam(form1_hour, data=df_train)
# summary(mod1_hour)
CV_gam(form1_hour, df_train, no_folds = 8)

##    score no_folds     exec_time
## 1 2.2825        8 1.913218 secs
##                                                                                                           formula
## 1 rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(hour, by = weekendsIndicator)

4 Modèles

On présente maintenant les différents types de modèle construits. Pour chacun, on construit un modèle avec et sans arêtes voisines, afin d’observer si ces dernières fournissent de l’information que le modèle convertit en une meilleure prédiction.

4.1 Arbres

Le premier modèle construit est un arbre aléatoire à l’aide du package rpart. On optimise le paramètre cp (complexity parameter) sur 2 grilles : la première de 0.0 à 0.1 et la seconde, plus fine, de 0.0 à 0.01.

Figure 4.1: Recherche du paramètre cp optimal

Le paramètre cp est le paramètre le plus important à optimiser. Dans la construction de l’arbre, si le gain de performance après une découpe n’est pas meilleur d’un facteur de cp, alors la découpe n’est pas réalisée. Dans notre cas, la valeur optimale est cp = 0, c’est-à-dire que l’on accepte toutes les découpes et que notre arbre est un arbre profond. On n’optimise pas les autres paramètres que l’on maintient à leur valeur par défaut car cp contrôle déjà la profondeur de l’arbre et on ne cherche pas à optimiser au maximum les arbres. On obtient les erreurs de test suivantes:

Table 4.1: Performances des arbres aléatoires
Arbre avec voisins	Arbre sans voisins
3.783	3.582

Premièrement, on note une performance significativement meilleure que les modèles naïfs. Ces arbres, bien que très simples, sont donc pertinents. Ensuite, on remarque que si on ajoute les arêtes voisines, le score est légèrement moins bon.

4.2 Forêts aléatoires

On construit maintenant des forêts aléatoires, une méthode d’ensemble à base d’arbres de décision. On utilise le package ranger pour les construire et le package vip pour visualiser l’importance des variables. On optimise 2 paramètres :

min.node.size (minimal node size) : l’algorithme ne découpe pas une feuille dont le nombre d’éléments est inférieur ou égal. Ce paramètre contrôle la taille des arbres.

Figure 4.2: Recherche du paramètre min.node.size optimal

On remarque que l’erreur décroît jusque la valeur 0 (en omettant la valeur au point 9). Pour ne pas que les arbres soient trop profonds, on garde la valeur par défaut \(\textbf{min.node.size} = 5\).

mtry (randomly selected predictors) : nombre de variables, tirées au hasard, que l’algorithme peut découper. Ce paramètre contrôle la corrélation entre les arbres.

Figure 4.3: Recherche du paramètre mtry optimal

Exceptionnellement, on affiche ici la recherche du paramètre optimal par la méthode timeslice. On remarque que l’allure globale des courbes sont identiques mais que le coude est plus prononcé avec timeslice. On choisit de poser \(\textbf{mtry} = 10\).

Le dernier paramètre à déterminer est ntree, le nombre d’arbres dans la forêts. Ce n’est pas un paramètre à optimiser comme les 2 précédents car le surapprentissage dans les forêts aléatoires n’est pas lié au nombre d’arbres. On optimise ntree afin de déterminer à partir de quel nombre d’arbre on gaspille du temps d’exécution pour un faible gain de performance.

Figure 4.4: Recherche du paramètre ntree ‘’optimal’’

On voit que le temps d’exécution est linéaire selon le nombre d’arbres et que la courbe de l’erreur marque un coude à partir de \(150\), valeur que l’on choisit pour **ntree*. On crée les forêts aléatoires et on obtient

Table 4.2: Performances des forêts aléatoires
Forêt aléatoire avec voisins	Forêt aléatoire sans voisins
3.092	3.129

Les performances des forêts aléatoires sont meilleures que celles des arbres aléatoires. Contrairement à ces derniers, les performances des forêts sont meilleures lorsque les arêtes voisines sont disponibles. De plus, on voit sur les scores d’importance qu’elles sont déterminantes dans les forêts. Par exemple, pour l’arête pont amont - pont austerlitz,

Figure 4.5: Scores d’importance dans l’ordre décroissante (15 premières arêtes) pour l’arête ‘’pont amont - pont austerlitz’’

Dans la suite, on souhaite construire des modèles nécessitant une sélection de paramètres. On la réalise grâce notamment grâce à ces scores d’importance. On souhaite également savoir quelles prédictions sont meilleures selon qu’on ajoute ou non les voisins. On représente sur le graphique suivant les différences entre les erreurs RMSE avec voisins et sans voisins. Le trait vertical noir point-tillé est la séparation entre Paris et le périphérique. Le trait horizontal rouge point-tillé est un seuil arbitraire (fixé à 0.2) à partir duquel on considère qu’ajouter les voisins améliore la performance de la forêt. Les arêtes respectant le seuil sont représentées par une croix rouge, sinon par une croix noire.

Seules 12 forêts sur 69 ont une performance ‘’significativement’’ meilleure lorsqu’on ajoute les voisins. De plus, pour les arêtes du périphérique, on remarque que les différences sont plus importantes que pour Paris, qu’inclure les voisins est néfastes (à 2 exceptions près) et qu’en général, la circulation sur ces arêtes est moins bien prédite. En effet, les scores pour Paris sont 2.747 (avec voisins) et 2.898 (sans voisins), et les scores du périphérique sont 4.069 (avec voisins) et 3.787 (sans voisins).